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^This  is  the  final  report  describing  a  throe  year  theoretical  investigation 
of  the  dynamics  of  fluid  spin-up  in  a  partially-filled  cylindrical  cavity.  The 
effort  was  divided  into  two, phases.  TTu*  first  phaso  consisted  of  developing  an 
extension  of  the  earlier  analyses  of  Nedem eyer„  and  of  Goller  and  Ranov,  to  those 
cases  where  the  liquid  free  surface  intersects  one  or  both  endwalls.  The  sispli- 
fying  assumptions  of  a  columnar  flow  and  the  quasi-steady  treatment  of  the  Ekman 
layer  pumping  of  the  secondary  flow  are  retained.  Earlier  estimates  of  the  Ekman 
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layer  pumping  are  modified  heuristically  for  situations  where  the  layer (s)  no 
longer  covers  the  entire  wall.  Also,  due  to  the  very  steep  free  surface  contour 
in  the  latter  stages  of  spin-up,  it  was  f?und  advantageous  to  develop  the  free 
surface  equations  in  an  axial,  rather  than  radial,  coordinate  frame.  This  model 
and  the  resulting  computer  code  are  the  subject  of  an  earlier  interim  report  and 
a  subsequent  paper  submitted  for  journal  publication.  Accordingly,  only  a  summary 
of  this  phase  and  its  results  is  presented  here. 

Ko  experimental  or  numerical  data  were  available  against  which  to  compare 
the  simplified  model's  predictions  for  cases  where  the  free  surface  intersects 
one  or  both  endwalls.  Accordingly,  in  the  second  phase  of  the  program  a  more 
refined  numerical  model  was  developed,  based  on  the  full  nonlinear  Navier-Stokes 
equations.  This  is  a  finite  difference  code  in  which  the  primitive  variable 
equations  are  solved  in'  conservative  form.  The  development  of  the  governing 
equations,  and  the  semi-implicit  predictor-corrector  scheme  used  to  solve  them, 
are  presented  in  some  detail^  Unfortunately,  no  converged  solutions  have  been 
obtained  as  yet;  the  possible^ reasons  for  this  are  discussed. 
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1.  INTRODUCTION 


The  goal  of  the  theoretical  program  described  here  was  to  obtain  a  clearer 
understanding  of  the  fluid  dynamics  of  liquid-filled  shells  for  those  cases  where  the 
fluid  only  partially  fills  its  cavity.  Problems  can  arise  during  the  flight  of  such 
projectiles  when  there  is  a  strong  coupling  between  the  perturbed  motion  of  the  rapidly 
spinning  shell  casing  a  d  the  attendant  wave  motion  excited  in  the  fluid.  In  addition, 
unlike  a  solid  filler  which  always  rotates  with  the  same  angular  velocity  as  the  casing, 
the  liquid  filler  will  rotate  at  a  reduced  angular  velocity  over  much  of  the  trajectory. 
The  result  is  a  net  transfer  of  angular  momentum  from  casing  to  fluid,  which  in  itself 
can  be  destabilizing.1 2 3 4 5 6 7 


The  model  problem  most  often  studied  by  investigators  is  the  axisymmetric  spin- 
up  of  fluid  in  a  right  circular  cylinder  of  radius  R  (see  Fig.  1).  We  seek  to  represent 
the  evolution  with  time  of  the  fluid  motion,  particularly  the  azimuthal  velocity,  as  it 
asymptotically  approaches  solid-body  rotation.  Such  a  calculation  is  necessary  as  input 
for  the  prediction  of  the  spin  decay  of  liquid-filled  projectiles.?*^  It  can  also  serve  as 
the  base  flow  about  which  to  carry  out  a  perturbation  analysis  of  the  eigenf requencies 
of  the  fluid.  ^  The'  latter  are  of  critical  importance  to  the  question  of  the  stability 
of  the  fluid-shell  system.1 


A  review  of  previous  work  on  this  problem  has  been  given  in  Refs.  6  and  7  and 
will  not  be  repeated  here.  The  principal  conclusion  drawn  from  these  studies  is  that, 


1.  Engineering  Design  Handbook.  Liquid-Filled  Projectile  Design,  AMC  Pamphlet 
No.  706-165,  U.S.  Army  Materiel  Command,  Washington,  D.C.,  April  1969. 

2.  Kitchens,  C.W.,  Jr.,  Gerber,  N.  and  Sedney,  R.,  "Spin  decay  of  Liquid  Filled 
Projectiles,"  J.  of  Spacecraft,  Vol.  15,  No.  6,  348-354,  December  1978. 

3.  Kitchens,  C.W.,  Jr;  and  Gerber,  N.,  Prediction  of  Spin  Decay  of  Liquid  Filled 
Projectiles,  BRL  Report  1996,  Aberdeen  Proving  Ground,  Md.,  July  1977. 

4.  •  Kitchens,  C.W.,  Jr.,  Gerber,  N.,  and  Sedney,  R.,  Oscillations  of  a  Liquid  iri  a 

Rotating  Cylinder:  Part  I.  Solid-Body  Rotation.  ARBRL-TR-02081,  Aberdeen 
Proving  Ground,  Md.,  June  1978. 

5.  Sedney,  R.  and  Gerber,  N.,Oscillations  of  a  Liquid  in  a  Rotating  Cylinder:  Part 
II.  Spin-Up,  ARBRL-TR-02489,  Aberdeen  Proving  Ground,  Md.,  May  1983. 

6.  Homicz,  G.F.,  Numerical  Model  for  Fluid  Spin-Up  in  a  Partially-Filled  Cylinder. 
Calspan  Report  No.  6856-A-l,  May  1982,  prepared  for  U.5.  Army  Research  Office 
as  interim  report  on  Contract  DAAG29-81-C-0G07. 

7.  Homicz,  G.  F.  and  Gerber,  N.,  "Numerical  Model  for  Fluid  Spin-Up  from  Rest  in 
a  PartiaTy-Filied  Cylinder,"  submitted  to  ASME  Transactions,  J.  of  Fluids 
Engineering; 
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Figure  1  CYLINDRICAL  GEOMETRY  AND  POSSIBLE  FLUID  CONFIGURATIONS  DURING 
SPIN-UP 


at  large  Reynolds  number,  the  dominant  spin-up  mechanism  is  a  secondary  flow  driven 
by  the  viscous  endwall  Ekman  layers.  This  secondary  flow  redistributes  the  angular 
momentum  acquired  in  the  Ekman  layers  over  the  remaining  interior  "core"  flow,  which 
is  relatively  inviscid.  The  characteristic  spin-up  time,  ts,  has  been  found  to  scale  as 
Of/ft)  (see  Nomenclature  for  a  list  of  symbol  definitions).^  For  large  Re,  this 

is  much  shorter  than  the  time  one  would  expect  based  so'ely  on  viscous  diffusion  from 
the  sidewall,  to,  which  scales  as  Re  A*'  . 

The  great  majority  of  previous  investigations  dealt  with  the  case  where  the  fluid 
cavity  is  completely  filled;  in  practice,  however,  the  shells  are  often  only  partially 
filled.  Field  data^  show  that  the  probability  of  a  projectile  experiencing  an  erratic 
flight  is  a  strong  function  of  its  fili  ratio,  i.e.,  the  percentage  of  the  cavity  volume 
occupied  by  fluid.  Thus  there  is  an  important  need  to  understand  the  influence  of  the 
interior  free  surface  on  the  spin-up  process.  Gerber  10  has  shown  that  the  final  form 
of  this  sirface,  after  solid-body  rotation  is  achieved,  is  a  parabola  whose  shape  is 
determined  by  the  fill  ratio  and  the  Froude  number,  F  =  Cfl.R)^/gH.  Goller  and  Ranov^ 
studied  how  the  surface  contour  evolves  in  time  and  its  influence  on  spin-up.  They 
employed  a  simplified  model  which  extended  the  partly  phenomenological  analysis  of 
Wedemeyer  for  the  filled  case.  Their  analysis  indicated  that  the  free  surface  motion 
acts  to  retard  the  spin- up  process.  Numerical  results  for  both  the  azimuthal  velocity 
field  and  the  surface  contour  were  obtained  as  functions  of  time,  and  the  latter  exhibited 
good  agreement  with  their  experimental  data.  The  problem  with  applying  their  analysis 
to  actual  configurations  is  that  Goller  and  Ranov  assumed  the  free  surface  did  not 
intersect  either  of  the  endwalls.  This  is  a  rathe.'  restrictive  assumption,  as  shown  by 
the  analysis  of  Gerber. 

Accordingly,  the  goal  of  the  first  phase  of  the  present  investigation  was  to 
extend  the  simplified  Wedemeyer -Goller -Ranov  model  to  situations  where  the  free 

8.  Greenspan,  H.P.,  The  Theory  of  Rotating  Fluids,  Cambridge  University  Press, 
London,  1968. 

9.  Mark,  A.,  Measurements  of  Angular  Momentum  Transfer  in  Liquid-Filled 
Projectiles.  ARBRL-TR-2029,  Aberdeen  Moving  Ground,  Md.,  November  1977. 

10.  Gerber,  NM  "Properties  of  Rigidly  Rotating  Liquids  in  Closed  Partially  Filled 
Cylinders,"  ASME  Transactions,  3.  of  Applied  Mechanics,  Voi.  97,  734-739,  1975. 

11.  Goller,  H.  and  Ranov,  TM  "Unsteady  Rotating  Flow  in  a  Cylinder  with  a  Free 
Surface,”  ASME  Transactions,  3.  of  Basic  Engineering,  Vol.  90D,  No.  4,  445-454, 
December  1968. 

12.  Wedemeyer,  E.H.,  "The  Unsteady  Flow  Within  a  Spinning  Cylinder,"  3.  of  Fluid 
Mecharics,  Vol.  20  Pt.  3,  383-399,  1964. 


surface  may  intersect  one  or  both  endwalls.  This  portion  of  the  work  is  discussed  in 
Section  2.  Unfortunately,  there  appear  to  be  no  data  for  such  flows,  either  experimental 
or  numerical,  gainst  which  to  compare  the  simplified  model's  predictions.  So  in  the 
second  phase  of  the  investigation  we  undertook  the  development  of  a  more  refined 
model,  based  on  the  full  Navier-Stokes  equations,  which  removes  the  principal 
approximations  of  the  earlier  analysis.  The  form  of  the  equations  and  the  finite 
difference  algorithm  used  to  solve  them  are  described  in  Section  3.  Section  4  summarizes 
the  conclusions  drawn  from  the  investigation  and  suggests  where  further  work  is  needed. 


2.  SIMPLIFIED  MODEL 


This  model  has  been  the  subject  of  an  earlier  interim  report.^  Since  its 
preparation,  a  few  improvements  have  been  marie,  arid  the  revised  model  and  its  results 
are  the  subject  of  a  manuscript  submitted  for  journal  publication.?  Hence  only  a  brief 
summary  of  the  principal  approximations  and  equations  will  be  given  here. 

2.1  Principal  Assumptions  and  Equations 

We  wish  to  predict  the  spin-up  of  fluid  in  a  cylinder  of  radius  R  and  height  H 
(see  Fig.  1)  which  at  t  =  0  impulsively  begins  to  spin  with  constant  angular  velocity, 
A,  about  its  axis.  The  fluid  is  incompressible  and  characterized  by  its  density,  p  , 
and  kinematic  viscosity,  V  .  Its  initial  level  when  at  rest  is  denoted  by  L.  The  velocity 
components  in  the  (r,  •,  z)  coordinates  are  (u,v,w)  respectively.  Time  is  here  normalized 
byjT,"l,  velocities  by  AR,  and  the  radial  and  axial  dimensions  by  R  and  H,  respectively. 

As  noted  earlier,  the  present  model  builds  on  the  previous  investigations  of 
Wedemeyer^,  and  Goller  and  Ranov^,  In  the  fully  filled  case,  Wedemeyer  argued  that 
in  the  limit  of  large  Re  the  flow  could  be  conceptually  divided  into  two  regions:  the 
viscous  flow  in  the  endwall  Ekman  layers,  and  the  so-called  "core"  flow  away  from  the 
endwalls  where  viscous  effects  are  much  less  important.  In  this  limit,  the  u  and  w 
velocity  components  are  much  smaller  than  the  azimuthal  component,  v,  so  that  most 
terms  in  the  Navier -Stokes  equations  involving  u  and  w  are  neglected.  This  leads  to 
a  simplified  model  in  which  u,  v  and  the  pressure  p  are  functions  of  only  r  and  t,  the 
so-called  "columnar"  flaw  approximation.  The  momentum  equation  governing  v  (r,t)  in 


the  core 

flow '  reduces 

to, 
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The  free  surface  contour  is  denoted  by  Zps  t>.  Following  Goller  and  Ranov,U  its 
evolution  is  governed  by  an  equation  which  balances  the  centipetal  and  gravitational 
force  components  tangential  to  the  surfaces 

,  p  £_* 

dr  r  (3) 

Equation  (1)  is  solved  using  a  Crank-Nicholson  type  implicit  finite-difference 
scheme  on  a  uniform  radial  grid.  Equation  (3)  is  then  integrated  numerically,  using  a 
Simpson's  rule  quadrature,  to  update  the  surface  contour.  Note  that  Eq.  (1)  is  still 
nonlinear  owing  to  the  presence  in  the  convective  term  of  u,  which  is  itself  a  function 
of  v.  Hence  an  iterative  process  is  required  at  each  time  step  to  get  a  consistent 
solution;  the  details  can  be  found  in  Ref.  6. 

The  relation  between  u  and  v  is  the  missing  link  needed  to  close  the  above 
system.  For  the  f»Ued  case,  Wedemeyer^  derived  such  a  relation  by  arguing  that 
whatever  outward  radial  mass  flux  is  generated  in  the  endwail  Ekman  layers  must  be 
balanced  at  each  radial  station  by  an  equal  but  opposite  inward  flux  in  the  core  flow. 
The  flux  in  the  Ekman  layers  is  known  to  develop  on  a  time  scale  very  much  faster 
than  the  spin-up  time.  Hence  its  magnitude  can  be  predicted  at  each  instant  by 
assuming  that  it  responds  in  a  quasi-steady  manner  to  changes  in  the  azimuthal  velocity 
of  the  core  flow  above  it.  Since  the  columnar  nature  of  the  flow  dictates  that  the 
compensating  radially  inward  flux  in  the  core  flow  be  spread  uniformly  in  the  axial 
direction,  this  gives  the  needed  relation  between  u  and  v.  It  is  this  inward  flux  that 
is  responsible  for  distributing  the  angular  momentum  acquired  in  the  endwail  layers 
over  the  remaining  fluid. 

The  preceding  argument  gives  the  contribution  to  u.in  the  core  flow  caused  by 
Ekman  layer  pumping.  Goller  and  Ranov^  reasoned  that,  in  the  partially  filled  case, 
this  would  be  opposed  by  an  outward  flux  in  the  core  owing  jo  the  fact  that  fluid  is 
being  flung  out  toward  the  sidewall  and  away  from  the  centerline.  Specifically,  the 
time  rate  of  change  of  fluid  volume  between  the  centerline  and  the  cylindrical  surface 
at  any  radius  r  must  equal  the  flux  of  fluid  across  the  surface.  Again  the  columnar 
approximation  requires  that  this  be  spread  uniformly  in  the  axial  direction.  The  result 
is  an  expression  for  the  positive  contribution  to  u  created  by  the  free  surface  motion. 
To  the  extent  that  it  opposes  the  negative  contribution  driven  by  the  Ekman  layers, 
the  spin-up  process  is  retarded. 
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Tbe  specific  form  of  the  equations  for  u  (v)  may  be  found  in  Ref.  6  and  7.  For 
flows  in  which  the  free  surface  does  not  intersect  either  endwall,  the  model  is  essentially 
the  same  as  Goller  and  Ranov's.  The  principal  difference  is  the  use  here  of  an  implicit, 
as  opposed  to  their  explicit,  time-marching  algorithm;  this  was  done  to  allow  the  use 
of  increased  time  steps,  and  hence  shorter  running  times. 


The  present  model  departs  from  that  of  Goller  and  Ranov  in  cases  where  one 
or  both  of  the  endwalls  are  intersected  by  the  free  surface.  Of  particular  interest  is 
the  question  of  how  best  to  describe  the  surface  contour  after  both  endwalls  are 
intersected.  Now  both  limits  of  integration  in  applying  Eq.  (3)  (r0  and  rj^  in  Fig.  1) 
become  functions  of  time  and  the  surface  contour  is  very  steep;  this  complicates  the 
accurate  integration  of  Eq.  (3)  with  respect  to  r.  Bo.h  difficulties  can  be  overcome 
by  switching  to  a  description  of  the  sirface  as  Rpg  (z,  t),  and  reinterpreting  Eq.  (3) 
with  z,  rather  than  r,  as  the  independent  variable  of  integration.  With  this  frame  of 
reference,  the  contour  has  a  shallow  slope  and  lies  between  constant  limits  of  integration. 
The  price  one  pays  is  that  an  integral  equation  must  be  inverted  to  obtain  Rp$  (z,  t), 
but  this  turns  out  to  be  amenable  to  a  simple  relaxation  solution  which  usually  converges 
in  a  few  iterations.  For  details,  the  reader  is  referred  to  Refs.  6  and  7. 


Another  distinguishing  feature  of  the  present  model  is  that,  once  the  free  surface 
has  intesected  an  endwall,  the  Ekman  layer  there  no  longer  completely  wets  the  surface. 
To  the  author's  knowledge,  there  have  been  no  published  investigations  into  the  structure 
of  the  layer  which  forms  in  such  a  case.  Accordingly,  heuristic  reasoning  was  used  to 
suggest  ample  analytical  modifications  to  the  Wedemeyer  expressions  which  reflect  the 
expected  reduced  mass  flux,  and  still  remain  within  the  columnar  flow  approximation. 
These  expressions  were  quoted  in  Ref.  6  as  Eq.  (19)  for  Utf;.  and  Eq.  (25)  for  Ujjel, 
which  represent  t  net  contributions  to  u  in  the  core  flow  from  the  top  and  bottom  Ekman 
layers,  respectively. 


Z2  ModificatiOi 


ns  To  Original  Model 


Subsequent 
to  better  reflect 
heuristic  modeling 
described  above, 
the  top  endwall  \ 


to  the  preparation  of  Ref.  6,  modifications  were  made  to  the  analysis 
the  physical  phenomena.  One  of  these  changes  has  to  do  with  the 
of  the  Ekman  layer  pumping  from  a  partially  wetted  surface,  as  just 
The  expression  now  used  for  the  contribution  to  u  from  the  layer  on 


o 


u. 


tel 


0  o  *  r  4  rH 

a  _  =  -  (R/H)rKr  (i  -  rH/r/r{(u,) 


rH  6  r*  i 


(4a) 

(4b) 


which  replaces  Eq.  (19)  of  Ref.  6.  Here  r{-|(t)  is  the  radius  at  which  the  top  wall  is 
intersected  by  the  free  surface,  and  Kj  and  Iff  are  adjustable  constants.  The  quantity 
f  (<*>)  represents  the  dimensionless  mass  flux  integral  appearing  in  the  Wedemeyer  12 
model,  as  a  function  of  the  local  angular  velocity,  6)  =  v/r  (see  Eq.  (10)  of  Ref.  6). 

This  expression  satisfies  the  following  constraints:  u  remains  continuous  across 
r  s  rfj;  the  contribution  from  uj£j_  vanishes  at  both  r  =  r^j  and  1,  as  it  should;  the 
influence  of  the  top  endwall  Ekman  layer  will  increase  as  rj-j  decreases,  which  is 
intuitively  satisfying;  and  the  columnar  nature  of  the  flow  is  maintained.  It  also  reduces 
to  the  form  of  the  original  Wedemeyer  model  in  the  fully  filled  limit,  =  0;  this 
last  point  was  not  tr- «  of  Eq.  (19)  in  Ref.  6.  Practically  speaking,  however,  this 
change  had  a  negligible  influence  on  the  numerical  results  for  the  conditions  studied  here. 


An  expression  analogous  to  Eq.  (4)  above  was  used  for  urft  ,  and  satisfies  similar 
criteria  (Eq.  (25)  of  Ref.  6).  Certainly  other  analytical  forms  can  be  found  which 
would  satisfy  these  constraints.  But  in  the  absence  of  any  experimental  data  to  serve  as 
a  guide,  it  was  decided  to  employ  as  simple  a  form  as  possible. 

Another  modification  to  the  original  model  concerns  the  appropriate  boundary 
condition  to  apply  to  Eq.  (1)  once  a  portion  of  the  bottom  wall  is  exposed.  After  this 
occurs,  one  cannot  expect  Eq.  (2b)  to  apply,  as  there  is  no  longer  any  fluid  at  the  axis. 
At  the  contact  line,  r  =  rQ  in  Fig.  1,  one  is  faced  with  the  apparently  conflicting 

requirements  that  there  be  no  slip  at  the  solid  surface,  but  that  the  contact  line 

nevertheless,  be  allowed  to  move.  The  resolution  of  this  paradox  is  still  an  active  area 
of  research,  see  e.g..  Ref.  13,  but  beyond  the  scope  of  the  present  investigation.  In 
our  simplified  model  we  initially  chose  to  use  a  no-slip  condition,  via.  v(r0)  *  r0  (Ref. 
6).  This  produced  a  discontinuous  change  in  angular  velocity  near  the  origin  as  the 

free  surface  passed  through  this  point.  An  argument  cat  also  be  made  that  one  has 

no  right  to  impose  a  no-slip  condition  on  the  (relatively)  Invlscid  core  flow  described 


13.  Pismen,  L.M.  and  Nir,  A.,  "Motion  of  a  Contact  Line,"  Physics  of  Fluids,  Vol. 
25,  No.  1,  3-7,  January  1982. 
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by  Eq.  (1),  where  the  terms  needed  to  define  the  Ekman  layer  structure  have  already 
been  eliminated. 


•Accordingly,  an  alternative  boundary  condition  on  v(rQ)  is  now  used,  one  based 
on  the  kinematic  condition  that  a  particle  on  the  free  surface  must  remain  there.  It 
is  most  easily  derived  by  taking  the  integral  of  Eq.  (3)  from  r0  to  some  general  positon 
r,  differentiating  the  result  with  respect  to  time,  and  then  evaluating  it  at  r  =  r0(t). 
The  boundary  condition  one  obtains  is: 


vcv  - 

.  L  Fdr./dt  J 


»/2 


which  replaces  Eq.  (20)  of  Ref.  6. 


(5) 


2.3  Results 


After  changing  the  analysis  as  described  in  Section  2.2,  the  numerical  results 
presented  in  Ref.  6  were  rerun  with  the  modified  code.  These  include  cases  for  a  range 
of  Froude  number,  Reynolds  number,  and  fill  ratio.  The  revised  .esults  are  presented 
in  Ref.  7.  A  representative  sample  case  is  displayed  here  in  Fig.  2.  This  is  referred  to 
as  Case  3  in  Ref.  6  and  7  and  corresponds  to  Re  =  1.172  x  105,  F  =  3.5,  H/R  *  3.0,  and 
L/H  =  0.6.  The  azimuthal  velocity  profiles  vs.  radius  are  shown  in  Fig.  2a.  The  first 
nine  profiles  are  for  dimensionless  times  t  =  400  to  3600  in  increments  of  400.  The 
last  profile,  for  ts  =  4745,  integrates  to  a  total  angular  momentum  about  the  axis 
greater  than  99%  of  its  final  value,  at  which  point  the  calculation  was  stopped.  The 
fluid  then  is  for  all  practical  purposes  in  solid-body  rotation.  The  last  five  profiles  in 
Fig.  2a  end.  abruptly  on  the  left  end  at  the  point  r  =  r^Ct),  where  the  free  surface  has 
intersected  the  bottom  wall.  Since  there  is  no  fluid  to  the  left  of  this  point,  v  is 
undefined  there. 

This  is  more  clearly  seen  in  Fig.  2b,  which  shows  the  free  surface  profiles  for 
the  same  times.  The  asymptotic  nature  of  the  approach  to  solid-body  rotation  is 
evident  in  both  figures  from  die  closer  spacing  of  the  profiles  later  in  the  calculation. 
Note  also  that  the  surface  contours  become  quite  steep,  and  thus  warrant  the  extra 
effort  spent  in  switching  from  a  radial  to  an  axial  grid  for  the  surface  integrals  when 
both  endwalls  are  intersected. 


A  semi-logarithmic  plot  of  the  dimensionless  angular  momentum  deficit, 
l-L2(tVL2(°°)»  vs.  Re-l/2t  is  presented  for  Case  3  in  Fig.  3.  The  reason  for  scaling 
the  time  axis  in  this  way  is  because  linearized  analysis  indicates  that  the  dimensionless 
spin-up  time  scales  as  Re^/2  (Ref.  8).  For  comparison  purposes,  the  results  from  Cases 
4  and  5  (Ref.  7)  are  also  presented  on  this  plot;  these  cases  differ  from  Case  3  only  in 
the  values  chosen  for  Re,  which  are  indicated  on  the  figure.  For  each  case  two  flagged 
symbols  are  shown.  The  first  of  these  represents  the  instant  at  which  the  free  surface 
first  hits  the  top  wall;  the  second,  the  instant  at  which  it  intersects  the  bottom  wall. 

The  most  noteworthy  feature  of  Fig.  3  is  that  up  to  the  instant  where  the  bottom 
wall  is  intercepted,  the  data  in  each  case  lie  on  a  straight  line.  This  is  consistent  with 
the  finding  of  linearized  theory  for  the  fully  filled  case  that  the  time  dependence  of 
the  spin-up  process  is  a  simple  exponent  .al.  What  is  surprising  is  the  degree  to  which 
this  remains  true  here,  where  nonlinearities  and  a  free  surface  are  both  present.  After 
the  bottom  endwail  is  intercepted,  however,  the  numerical  data  in  Fig.  3  follow  a  line 
with  a  reduced  slope,  indicating  a  slower  rate  of  spin-up. 

The  fact  that  the  surface  hitting  the  top  wail  apparently  has  little  effect  on 
the  rate  of  spin-up  .can  be  explained  by  noting  that  the  radial  extent  of  the  fluid  at  the 
top  wall  in  Fig.  2b  is  relatively  thin;  moreover,  the  fluid  in  this  region  has  already 
accpjired  a  significant  fraction  of  its  final  angular  velocity.  Hence  its  contribution  to 
the  Ekman  pumping  remains  small  compared  with  that  from  the  bottom  layer,  which 
still  completely  covers  the  wall.  Once  the  bottom  wall  is  intercepted  in  the  final 
stage,  however,  the  strength  of  the  bottom  Ekman  layer  pumping  is  progressively  reduced 
as  the  wall  is  exposed,  and  a  slower  spin-up  rate  results.  Thus,  it  is  the  bottom  Ekman 
layer  which  is  the  dominant  driving  force  behind  the  secondary  flow  responsible  for 
spin-up. 

As  noted  in  Section  2.2,  Cases  3-3  had  originally  been  run  (Ref.  6)  with  a  no¬ 
slip  condition  on  v(r0),  whereas  the  results  reported  here  and  in  Ref.  7  used  Eq.  (5). 
A  comparison  of  the  two  sets  of  results  shows  that  the  latter  lowers  the  velocity 
profiles  !n  the  immediate  vicinity  of  r0,  and  increases  the  spin-up  time  by  3%,  5%  and 
19%  for  Cases  3-3,  respectively.  Referring  again  to  Fig.  3,  the  only  qualitative  change 
of  any  note  was  that,  with  the  no-slip  condition,  Case  3  had  previously  exhibited  no 
kink  as  the  bottom  wall  was  intersected.  This  was  interpreted  as  meaning  the  Reynolds 
number  in  this  case  was  low  enough  for  viscous  diffusion  from  the  sidewall  to  dominate 
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over  Ekman  layer  pumping.  With  the  revised  results  all  that  can  be  said  is  that  the 
change  in  rate  of  spin-up  at  the  second  transition  decreases  with  decreasing  Re  in  Fig. 
3.  Finally,  it  is  interesting  to  note  that  the  angular  momentum  ratios  at  the  times 
at  which  transitions  occur  are  nearly  independent  of  Re. 

Additional  numerical  results  and  discussion  may  be  found  in  Refs.  6  and  7. 
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3.  NAVIER-STOKES  MODEL 


The  simplified  model  described  in  the  preceding  section  is  predicated  or.  the 
following  principal  approximations:  conceptual  division  of  the  flow  into  the  viscous 
endwall  Ekman  layers,  and  the  relatively  inviscid  "core11  flow;  the  assumption  of  columnar 
flow  in  the  core;  and  the  quasi-steady  treatment  of  the  Ekman  layer  pumping.  Each 
of  these  becomes  more  questionable  as  the  Reynolds  number  is  lowered.  Moreover, 
the  only  experimental  data  available  to  validate  the  simplified  model's  predictions  are 
those  in  Ref.  11,  which  pertain  only  to  situations  where  the  free  surface  intersects 
neither  of  the  endwalls. 

Hence  the  goal  of  the  second  phase  of  the  present  program,  described  below, 
was  to  develop  a  more  refined  numerical  model  based  on  the  full  Navier-Stokes  equations. 
Such  a  model  would  relax  all  the  above  assumptions,  and  serve  to  validate  over  what 
range  of  parameters  the  simplified  model  could  be  trusted. 

3-1  Governing  Equations 

The  flow  in  principle  is  governed  by  the  continuity  equation, 

V  •  V  *  o  (6) 


f”  H"  ^  ‘  R*'  7  *  7  x  *  (7) 

where  the  time  t  is  again  normalized  by  the  velocities  by iXR,  and  the  pressure  by 
p  (AR&  in  contrast  to  the  simplified  model,  however,  all  lengths  are  here  normalized 
by  the  cylinder  radius,  R.  In  dimensionless  coordinates  then,  the  cylinder  maps  to  the 
region  0  4  r  4  1,  0  4  z  4  H/R. 

The  convective  term  in  Eq.  (7)  has  been  written  in  conservative  form  using  the 
identity 

if  V  V  *  V  ■  (is  S’)  -  v  ( 7  ■  v ) 


and  the  Navier-Stokes  equations. 


d  ir 

+  7  •  (V  is)  -  -  V  -p  - 

o  X 


while  the  viscous  term  was  transformed  using 

7*  v  *  -Vx(Vxir)  +  V(V-v) 


The  second  term  in  each  identity  vanishes  as  required  by  Eq.  (6),  and  the  result 
is  Eq.  (7).  In  this  form  the  finite  difference  analog  of  the  convective  terms  conserves 
momentum.  The  reason  for  transforming  the  viscous  terms  in  this  way  will  become 
apparent  shortly. 

Equations  (6)  and  (7)  comprise  four  scalar  equations  for  the  unknowns  u,  v,  w  and 
p.  An  additional  relation  is  needed  to  predict  the  unknown  free  surface  contour.  As 
before,  we  describe  the  free  sirface  by  G  =  z  -  Zps  (r,  t)  =  0.  The  kinematic  condition 
that  a  fluid  particle  on  the  surface  must  remain  there  then  requires  that  ** 

2JL  »  +  ir  -  VG  -  0  W 

Dt  at 

at  the  surface,  G  =  0. 

When  written  out  in  cylindrical  coordiarttes,  Eqs.  (6)  -  (8)  become 
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where,  e.g.,  wps  represents  the  axial  velocity  component  at  a  point  on  the  free  surface. 
The  only  assumptions  that  have  been  made  in  writing  the  above  equations  are  that  the 
flow  remains  axisymmetric  and  laminar. 

To  uniquely  define  the  problem,  initial  and  boundary  conditions  for  Eqs.  (6)  -  (8) 
are  needed.  They  are  most  easily  discussed  below  in  connection  with  the  finite  difference 
form  of  the  equations. 

3.2  Coordinate  Transformations 

In  order  to  better  resolve  the  viscous  Ekman  layers  on  the  endwalls  and  the 
Stew  art  son  layer  at  the  sidewall,  a  coordinate  stretching  is  employed.  We  use  the 
same  transformation  successfully  used  by  Kitchens ^  for  the  filled  case.  The  radial 
coordinate  r  is  transformed  to  p  where 


'  '  *>  (&)/*-(&) 


The  inverse  transformation  is 


Kitchens,  C.W.,  Jr.,  "Navier-Stokes  Solutions  for  Spin-Up  in  a  Filled  Cylinder," 
AIAA  Voi.  18,  No.  8,  9 29-934,  August  1980. 
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Similarly,  the  axial  coordinate  is  transformed  according  to 


/C  +  if  z  '  1 

whose  inverse  is  given  oy, 


The  above  transformations  map  the  cylinder  into  the  unit  square  0  £  p  £  i,  0  £  'f  £  1  in 
the  (p,  f  )  computational  plane.  The  corsiants  b  and  c  (41)  are  input  parameters  used 
to  adjust  the  grid  spacing;  as  b  and  c  approach  unity  from  above,  progressively  more 
points  will  be  clustered  near  the  side-  and  endvalls  in  the  physical  (r,  z)  plane, 
respectively. 

Upon  application  of  the  above  transformations,  Eqs.  (6H8)  become,  in  the 
computational  plane: 
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where  subscripts  are  used  to  denote  differentiation. 
3.3  Computational  Grid 


The  above  equations  are  differenced  on  a  uniform  rectangular  grid  in  the 
computational  (9 ,  ^  )  plane.  Initially  an  attempt  was  made  to  use  a  conventional  grid, 
i.e.  one  in  which  the  velocity  components  and  pressure  are  all  defined  at  the  same 
grid  locations.  But  it  was  found  that  this  leads  to  an  inconsistency  between  the 
difference  forms  of  the  velocity  and  pressure  equations.  This  will  be  mentioned  again 
below  in  connection  with  the  form  of  those  equations. 

The  differencing  scheme  now  used  employs  a  so-called  "staggered"  or  interlocking 
grid,  which  was  first  proposed  by  Welch,  et  al.16  in  connection  with  the  Marker-and- 
Cell  (MAC)  method.  Although  the  explicit  MAC  method  is  not  used  here,  for 
incompressible  flows  the  staggered  grid  arrangement  still  has  clear  advantages.  This 
arrangement  is  shown  in  Fig.  4.  The  letters  }  and  k  are  used  as  indices  in  the  radial 
and  axial  coordinates,  respectively.  The  pressure  p  is  defined  at  the  center  of  each 
grid  cell,  pj,  |<,  i«e.  where  integral  values  of  j  and  k  intersect.  The  radial  and  azimuthal 
velocities  are  defined  at  the  center  of  the  cell's  vertical  faces,  at  (j  >  1/2,  *);  and 
the  axial  velocity  is  defined  at  the  center  of  the  horizontal  faces,  at  (j,  k  ♦  1/2). 


16.  Welch,  J.E.,  Harlow,  F.E.,  Shannon,  3.P.  and  Daly,  The  MAC  Method?  A 
Computing  Technique  for  Solving  Viscous,  Incompressible,  Transient  Fluid-Flow 
Problems  Involving  Free  Surfaces,  Los  Alamos  Scientific  Laboratory  Report  LA- 
3423,  March  1966. 


Vi 


The  grid  is  numbered  in  such  a  way  that  the  cylinder  axis  and  the  sidewall  lie 
along  j  =  3/2  and  j  =  NJ  -  1/2,  respectively;  NJ  is  the  chosen  number  of  radial  grid 
cells.  The  relation  between  p  and  j  is  thus  J)  =  (j  -  3/2)  AjB  ,  where  =  l/(NJ-2). 
Analogously,  the  bottom  and  top  endwalls  lie  along  k  =  3/2  and  k  =  NIC  -  1/2,  respectively, 
with  NK  the  number  of  axial  cells.  Thus  f  =  (k  -  3/2)  A  f  ,  where  =  l/(NK-2). 

3.4  Difference  Equations  For  Velocities 

The  algorithm  used  to  solve  the  equations  is  an  adaptation  of  the  semi-implicit 
Predictor-Corrector  Multiple-Iteration  (PCMI)  method.  This  method  was  first  proposed 
by  Rubin  and  Lin  17  for  solving  steady,  three-dimensional  boundary  layer  problems. 
Kitchens  used  it  to  study  the  spin-up  problem  in  a  filled  cylinder,  which  is  believed 
to  be  the  method's  first  application  to  the  full,  unsteady  Navier-Stokes  equations. 
Because  of  the  absence  of  a  free  surface,  Kitchens  was  able  to  more  easily  formulate 
his  equations  in  terms  of  the  stream  functioh-vortidty  variables.  The  present  study 
is  thus  believed  to  be  the  first  time  the  PCMI  scheme  has  been  applied  to  the  primitive 
variable  formulation  of  the  Navier -Stokes  equations,  with  a  free  surface  present. 


The  basic  idea  of  the  PCMI  scheme  is  rather  simple.  We  will  use  a  uniform 
time  step,  At,  and  a  superscript  i  to  denote  the  time,  level  such  that  t  *  iAt.  At 
the  beginning  of  the  step  to  time  (i  +  1),  all  flow  variables  at  level  i  will  be  known. 
The  first  "predictor"  part  of  the  calculation  consists  of  extrapolating  all  flow  variables 
to  level  (i  ♦  1)  using  information  from  the  three  previous  level* 
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i+t 


3y‘"+ 
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which  is  accurate  to  order  (At)3.  Here  Y  is  a  generic  symbol  for  u,  v  or  w,,  and 
spatial  indices  have  been  dropped  for  simplicity.  For  the  first  two  time  steps, 
extrapolations  of  order  At  and  (At)2  are  used,  respectively.  Equrtion  (14)  provides 
starting  values  for  each  of  the  variables.  This  is  followed  by  an  iterative  process  in 
which  solution  of  the  Navier-Stokes  equations  serves  to  "correct*  the  initial  guess. 


17.  Rubin,  S.  G.  and  Lin,  T.  Cn  "A  Numerical  Method  for  Three-Dimensional  Viscous 
Flow*  Application  to  the  Hypersonic  Leading  Edge,"  J.  of  Computational  Physics, 
Voi.  9,  339-364,  1972. 


The  time  derivative  in  Eq.  (12)  is  approximated  at  level  (i  +  1)  by  second-order 
backward  differences, 

-  (3Yi*'-4'fi  +  Yi")/2&t 

For  consistency,  the  spatial  derivatives  on  the  right  side  of  Eq.  (12)  must  also  be 
evaluated  at  (i+1),  hence  the  implicit  nature  of  the  scheme.  For  example,  let 

denote  the  +  ^  iterate  to  the  «tesire<l  *  The  viscous  30(1 

pressure  gradient  terms  are  then  approximated  by  second-order  centered  differences. 
Derivatives  in  the  ya  coordinate  are  treated  implicitly,  i.e.  in  terms  of  (n  +  1)  iterate 
values.  Derivatives  in  the  f  direction  use  these  values  as  they  become  available,  with 
(n)  iterate  values  used  where  they  are  not,  as  the  grid  is  swept  from  bottom  to  top. 
For  example,  the  viscous  terms  in  Eq.  (12a)  are  approximated  by 


i*t 
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The  nonlinear  convective  terms  in  the  equations  are  differenced  using  a  hybrid 
of  centered  and  upwind  schemes,  as  suggested  by  Hirt,  et  al.^  As  an  example,  the 
radial  convection  term  in  Eq.  (12a)  becomes 


Hirt,  C.W*  Nichols,  B.D.  and  Romero,  N.C.,  SOLA  -  A  Numerical  Solution 
Algorithm  for  Transient  Fluid  Flows,  Los  Alamos  Scientific  Laboratory  Rejjort. 
LA-5852,  April  1975. 
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where  ot  is  an  adjustable  input  parameter  between  zero  and  one.  For  aC  =  Q,  centered 
differencing  is  recovered,  while  oc  =  1  yields  full  upwind  or  "donor-cell'*  differencing.! 9 
Such  a  scheme  is  stable  and  conservative,  yet  reduces  the  formal  truncation  error  to 
first-order;  as  pointed  out  in  Refs.  19  and  20,  it  may  still  produce  results  comparable 
in  accuracy  to  those  of  a  second-order  centered  scheme.  Note  that  the  convective 
velocity  is  evaluated  from  the  previous  iterate,  so  that  the  algebraic  equations  are 
rendered  linear.  The  above  form  calls  for  values  of  u  and  w  which  are  not  defined  in 
the  grid  system  of  Fig.  4.  For  such  quantities,  simple  averages  are  taken,  e.g., 


This  scheme  leads  to  a  tridiagonal  system  of  equations  for  each  of  the  velocity 
components  along  a  given  k  row.  Because  of  their  length,  the  full  system  is  given  in. 
the  Appendix.  For  each  k,  the  equations  are  easily  solved,  using  a  standard  algorithm 
for  inverting  tridiagonal  systems.!9  The  initial  conditions  at  t  a  0  are  u  =  v  =  w  = 
0  everywhere.  The  grid  is  swept  by  first  correcting  the  values  adjacent  to  the  endwall, 
k  *  2,  and.  proceeding  one  row  at  a  time  up  to  the  free  surface. 


At  the  cylinder  axis,  the  following  boundary  conditions  are  imposed: 
IL. 


•*,*  '  V\,* 


(15) 


19.  Roache,  P.J.,  Computational  Fluid  Dynamics.  Hermosa  Publishers,  Albuquerque, 
NJWl.,  1976. 

20.  Spalding,  D.  B.,  "A  Novel  Finite  Difference  Formulation  for  Differential 
Expressions  Involving  both  First  and  Second  Derivatives,"  Interna tionai  Journal 
for  Numerical  Methods  in  Engineering,  Voi.  4,  551-559,  1972. 
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At  the  bottom  wall,  a  no-slip  condition  is  enforced: 
iX  •  , '  =  -U.  ■  .j.  - 

V.  +x  i  *  2  Cl  r  -j_  -zr  +x  ' 

y-t.  y*-,* 

(Sr  i  *  o 
i>  x. 

And  at  the  sidewall,  the  following  obtain: 


**J-  x,* 
fN3-±,Je 


*  -Q 


where  2L  denotes  the  cylinder's  angular  velocity  normalized  by  its  final  value.  For 
ideally  impulsive  spin-up,  Jt  =  1  for  t  >  0.  The  code  allows  for  a  "ramped"  imposition 
of  this  boundary  condition,  i-e.. 


SL  -  t/t # 
SI  '  -  1 


t  4  t. 
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where  typically  t{Mp  a  10.  Equations  (15)  -  (17)  are  accurate  to  second  order. 

Finally,  there  is  die  question  of  what  boundary  conditions  are  appropriate  at  the 
free  surface.  Let  KT  (j)  denote  the  topmost  cell  in  the  j-column  which  contains  fluid. 
Following  Hirt,  et  al.1*  if  KT  Q)  does  not  exceed  KT  (j-1)  the  azimuthal  and  radial 
velocities  at  (j.  ♦  1/2,  KT),  along  with  all  other  interior  ceil  values  to  their  right,  are 
corrected  via  the  Navier-Stokes  equations  as  described  above;  otherwise,  they  are  set 
equal  to  the  values  in  the  cell  immediately  below.  In  either  case,  the  axial  velocity  in 
the  surface  cells  is  chosen  so  as  to  satisfy  the  difference  approximation  to  Eq.  (11),  i.e., 
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This  specification  assures  that  the  divergence  in  cell  (j,  KT)  vanishes  (see  Eq.  (20)  below). 


The  above  boundary  conditions  apply  to  situations  in  which  neither  endwall  is 
intersected  by  the  surface,  and  will  require  some  modification  for  cases  in  which  this 
is  not  so. 


The  semi-implicit  PCMI  method  described  above  was  chosen  over  simpler  explicit 
schemes,  e.g.  the  MAC  method, 16  for  several  reasons.  Since  radial  variations  are 
treated  implicitly,  the  allowable  At  depends  only  on  the  axial,  and  not  die  radial,,  grid 
spacing  (Refs.  15,  17).  Axial  gradients  in  the  core  flow  tend  to  be  much  smaller  (recall 
the  columnar  flow  approximation  in  Section  2),  and  so  will  allow  larger  grid  spacings 
and  time  steps.  Indeed,  this  is  the  reason  for  treating  the  radial  and  axial  variations 
implicitly  and  explicitly,  respectively.  The  PCMI  method  has  also  been  proven  capable 
of  handling  the  strong  nonlinear  coupling  between  the  equations,  as  demonstrated  by 
its  success  for  the  filled  cylinder  problem.^  Finally,  the  iterative  structure  of  the 
scheme  provides  a  convenient  framework  for  correcting ,  the  free  surface  contour  at 
each  step.  The  overall  iterative  process  will  be  discussed  in  Section  3.7. 


3.5  Difference  Equation  For  Pressure 

The  time  derivative  of  p  does  not  appear  in  Eqs.  (11W13).  Hence  they  cannot 
be  used  to  march  the  pressure  field  forward  in  time  in  the  same  manner  as  the  velocity 
field.  Instead,  the  pressure  field  is  adjusted  at  each  time  step  so  as  to  enforce  the 
incompressibility  condition  expressed  by  Eq.  (11).  Let  D  =  V*  v  ,  and  Dj^  denote  its 
value  at  the  point  (j,  k)  in  Fig.  4.  Using  centered  differences. 
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The  equation  for  p  is  obtained  by  taking  the  same  linear  combination  of  Eqs.  (A-l)  and 
(A-3)  as  appears  in  Eq.  (20).  One  can  see  by  inspection  that  the  terms  involving  time 
derivatives  will  lead  to: 


2  A£ 


(21) 


which  is  a  second-order  accurate  approximation  to 
Those  terms  involving  p  combine  to  yield 


X  L" 

3t  'j.-k 
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which  is  a  second-order  approximation  to  the  Laplacian, 
d-p  \ 


*t>r 


centered  at  (j,  k).  The  terms  which  result  from  the  convection  terms  in  the  momentum 
equations  are  too  lengthy  to  quote  here  in  theii  original  form.  However,  they  are 
easily  recognized  as  the  difference  approximation  to  the  continuum  term,  7* 0?) at 
(j,k).  This  brings  up  another  important  point  which  differentiates  the  present  scheme 
from  the  MAC  method.  One  can  easily  show,  using  tensor  notation  and  the  Einstein 
summation  convention,  that 
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The  last  two  terms  vanish  for  the  continuum  equations,  as  per  Eq.  (6).  In  the  explicit 
MAC  method,  all  three  terms  are  essentially  included  in  the  pressure  equation  in 
difference  form.  But  because  they  get  evaluated  at  the  previous  time  step,  for 
which  OjA  so  is  already  satisfied  to  some  prescribed  accuracy,  the  contribution 
from  the  last  two  terms  on  the  right  is  made  vanishingly  small  -  as  it  should  be. 

In  the  present  implicit  formulation  these  terms  are  now  evaluated  at  the  new 
time.  However,  in  general  will  not  be  small  during  the  iteration  process  prior 

to  convergence,  and  so  the  presence  of  the  last  two  terms  on  th.  rij-rt  can  be 

£♦  i 

destabilizing.  For  this  reason,  and  to  be  consistent  with  setting  D ^  =  0  in  the 

time  derivative  (see  below),  we  omit  the  last  two  terms  on  the  right  from  the  pressure 
equation.  Thus  the  cylindrical  coordinate  representation  of  |^.  is  all  that  remains 
of  the  convective  terms.  The  full  difference  equation  for  p  is  quoted  in  the  Appendix 
as  Eq.  (A-4). 

By  making  the  above  term-by-term  identifications,  we  see  that  the  substitution 
of  Eqs.  (A-l)  and  (A- 3)  into  Eq.  (20)  yields  exactly  the  same  difference  equation  as 
would  be  obtained  by  first  taking  the  divergence  6f  Eq.  (7), 

«  -££  -  v-(v-vv) 

9t 

and  then  applying  the  difference  approximations.  At  first  reading,  this  statement  might 
appear  obvious.  But  in  fact,  as  discussed  by  Welch,  et  al.,16  the  staggered  grid 
arrangement  of  Fig.  4  appears  to  be  the  only  one  for  which  this  statement  is  true.  In 
particular,  it  is  not  true  for  a  conventional  grid  arrangement,  for  which  the  two 
approaches  produce  entirely  different  results.  It  was  because  of  just  this  inconsistency 
that  the  conventional  grid  was  abandoned  here  in  favor  of  the  staggered  grid  (cf . 
Section  3.3). 

Another  noteworthy  feature  of  Eq.  (A-4)  is  the  complete  absence  of  any  viscous 
terms,  which  turn  out  to  be  self-cancelling.  This  is  the  reason  for  having  expressed 
the  viscous  terms  in  the  form  used  in  Eq.  (7).  Since  the  divergence  of  the  curl  of 
any  vector  field  is  identically  zero,  a  considerable  simplification  of  the  pressure  equation 
results. 
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Equation  (A-4)  recognizes  that  while  D  is  supposed  to  vanish  throughout  the  flow 
at  all  time  steps,  because  we  can  never  solve  the  equations  exactly  some  small  but 
finite  divergence  will  be  present  in  the  solutions  at  times  i  and  (i-1).  However,  the 
pressure  field  is  now  to  be  determined  in  such  a  way  as  to  make  the  divergence  vanish 
at  time  0+1).  This  is  reflected  in  Eq.  (A-4)  by  setting  only  the  term  from 

expression  (21)  to  zero,  while  the  values  at  i  and  (i-1)  are  found  from  the  velocity  field 
at  those  times.  In  this  way,  any  growth  in  the  divergence  field  is  self-limiting,  and 
the  solution  remains  stable.^ 

Equation  (A-4)  represents  a  Poisson  equation  foe  the  pressure  with  the  source 
term  on  the  right  a  function  of  the  latest  iterate  to  the  velocities  at  (i+1).  It  is 
solved  here  using  a  standard  point  SOR  method,  traversing  the  grid  starting  at  the 
lower  left  corner  (j  =  2,  k  =  2),  passing  from  left  to  right  across  each  row,  and  sweeping 
from  the  bottom  row  up  to  the  free  surface.  The  initial  condition  is  just  the  hydrostatic 
pressure  distribution. 


the  starting  value  for  which  is  given  by 


The  boundary  conditions  at  the  axis  and  solid  boundaries  are  obtained  by 
specializing  the  appropriate  normal  momentum  difference  equation  to  that  boundary, 
and  employing  Eqs.  (15)  -  (17).  Hence  at  the  centerline, 
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At  the  bottom  wall, 


At  the  sidewall. 


where  again  JT.  has  been  used  rather  than  its  final  value  of  unity  to  allow  for  the 
ramped  imposition  of  the  sidewall  boundary  condition,  cf.  Eq.  (18).  The  above  represent 
Neumarav-type  boundary  conditions  on  the  solution;  they  are  incorporated  directly  into 
the  relaxation  equations  at  the  points  immediately  adjacent  to  the  boundary  as  described 
by  Roache19,  pp.  183-184. 

The  pressure  in  the  surface  cell  at  (j,  KT(j»  is  found  not  from  Eq.  (A-4),  but 
rather  from  the  boundary  condition  that  the  pressure  at  the  free  surface  be  a  constant, 
say  pps*  This  boundary  condition  neglects  the  effects  of  surface  tension,  which  are 
small  for  the  surface  radii  of  curvature  of  interest  here.**  Also  neglected  are  viscous 
stresses  at  the  surface;  their  influence  has  been  shown  by  experience  to  be  significant 
only  when  Re  10.21  Referring  to  Fig.  5,  we  denote  by  h  the  distance  the  free  surface 
lies  above  the  pressure  node  at  KT,  expressed  in  the  transformed  f  coordinate.  Then 
a  linear  interpolation  between  the  surface  and  the  node  at  KT-l  yields 


*7,»r  -  O,,*  + 


21.  Hirt,  C.W.  and  Shannon,  J.P.,  "Free-Surface  Stress  Conditions  for  Incompressible- 
Flow  Calculations,”  3.  of  Computational  Physics,  Voi.  2,  No.  4,  403-411,  1963. 
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Figura  5  PRESSURE  BOUNDARY  CONDITION  AT  THE  FREE  SURFACE 
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Similarly,  if  the  value  at  KT+1  should  be  needed, 

1°j.Kr+i  s  2’f°j,KT  "  (25b) 

Hence  only  in  interior  cells  is  the  pressure  adjusted  to  satisfy  Eq.  (20);  in  surface  cells 
this  function  is  served  by  the  velocity  boundary  conditions,  Eq.  (19),  with  Eq.  (25a)  used 
to  determine  the  pressure. 


3.6  Difference  Equation  For  Surface  Contour 

The  continuous  function  Zpg  ( j> ,  t)  in  Eq.  (13)  is  discretized  by  defining  it  at 
integral  values  of  j,  i.e.  Zi^.  denotes  Zfs  at  /  3  f*\  ****  *  s  1  At*  Second-order 
backward  differences  are  again  used  for  the  time  derivative.  The  convective  term  on 
the  right  is  approximated  using  a  hybrid  of  upwind  and  centered  differences  analogous 
to  that  used  in  the  momentum  equations.^  The  result  is 


!>»♦* 


-2 
(26) 
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where  is  an  input  parameter  analogous  to  *  ;  values  of  ^3  =  0  and  1  produce  purely 
centered  and  upwind  differencing,  respectively.  zR§j  represents  the  (n+1)  iterate  to 
the  contour  ,  and  w£+I  the  free  surface  axial  velocity  at  j>  ,  z  a'  Zps.  found 
by  linear  interpolation  of1  the  (n+1)  velocity  iterate.  Eq.  (26)  represents  a J linear, 
tridiagonal  system  for  the  unknown  Z^f,  and  so  can  be  inverted  using  the  same  efficient 
algorithm  as  for  the  momentum  equations.  The  initial  condition  is  Zpgj.a  L/R  ait  t  *  0 
for  all  j. 


Since  Eq.  (13)  is  of  first  order  in  the  radial  coordinate,  only  a  single  boundary 
condition  is  required.  But  the  values  of  Zpg  at  the  centerline  and  the  sidewall  are 
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not  known  a  priori.  So  instead,  a  global  boundary  condition  is  employed,  based  on  the 
fact  that  total  fluid  volume  is  conserved.  At  each  instant  this  requires,  in  dimensionless 
form, 


,  v-l 

]jk  ri  2'H  ° 


The  above .  summation  is  a  trapezoidal  rule  approximation  to  the  volume  integral,  in 
terms  of  the  transformed  variable  p  .  The  constant  V0  is  easily  determined  at  the 
start  when  all  the  Zpsj  equal  the  initial  fill  level.  To  see  how  this  condition  is 
enforced,  we  note  that  the  last  step  in  the  tridiagonal  solution  algorithm  19  for  the 
Zn+lpsj  **  the  application  of  the  following  linear  recursion  between  neighboring  values: 

Snj  -  Sf  W) 

Ej  and  Fj  represent  coefficients  derived  from  Eq.  (26)  which  depend  on  the  surface 
velocities  and  p  ,  but  not  on  the  Zn+ipgj: 

From  Eq.  (2S)  we  see  that  once  the  fluid  level  at  j  *  NJ  is  given  any  value  the 
remainder  of  the  contour  is  uniquely  determined.  So  also  is  the  summation  in  Eq.  (27), 
which  in  general  would  not  equal  V0.  However,  since  Eq.  (28)  is  linear  in  Zpg  one 
can  easily  substitute  it  into  Eq.  (27)  and  solve  analytically  for  the  value  at  j.  a  N3,  in 
terms  of  VQ  and  other  constants,  all  of  which  are  already  known.  In  this  way,  the 
needed  boundary  condition  on  the  surface  contour  is  determined  so  as  to  automatically 
satisfy  conservation  of  fluid  volume  without  the  need  for  additional  iteration. 

3.7  Overall  Computational  Cycle 


Here  we  outline  the  principal  steps  involved  in  advancing  the  solution  through 
one  time  cycle  using  the  results  of  the  preceding  sections. 


It  is  assumed  that  the  solution  at  time  levels  1,  i-l»  and  i-2  is  known, 
either  from  the  initial  conditions,  or  because  we  have  marched  the 
calculation  that  far. 


A  prediction  is  made  for  each  of  the  velocity  components,  the  pressure, 
and  the  surface  contour  using  the  extrapolation  in  Eq.  (14).  This  can  be 
viewed  as  the  zeroth  iterate  to  the  solution  at  i+1. 

The  velocity  field  is  corrected,  i.e.  the  (n+1)  st  iterate  is  obtained  from 
the  solution  of  Eqs.  (A-l)  to  (A-3),  as  described  in  Section  3.4. 

The  (n+l)st  iterate  to  the  velocity  field  is  used  to  interpolate  for  the 
radial  and  axial  velocity  components  at  the  surface,  and  then  Eq.  (26)  is 
solved  for  the  (n+l)st  iterate  to  Zi+1FSj* 

The  boundary  conditions  on  velocity  and  pressure  across  the  new  iterate 
to  the  free  surface  are  updated  using  Eqs.  (19)  and  (25). 

Equation  (A-4)  is  solved  using  point  SOR  to  obtain  the  (n+l)st  iterate  to 
the  pressure  field  at  time  (i+1).  This  actually  represents  an  inner  iteration 
nested  inside  the  outer  iteration  loop  for  the  velocity  field.  The  pressure 
solution  is  considered  to  have  converged  when 


where  Apj|<  is  the  change  in  pressure  at  (j,  k)  from  one  pressure  iteration 
to  the  next.  €  rp  and  €ap  are,  respectively,  the  relative  and  absolute 
convergence  criteria  applied  to  the  pressure  solution;  typical  values  used 
are  €rp  a  10“*  and  €ap  3  10-6. . . 

The  pressure  boundary  condition  across  the  free  surface  is  updated  using 
Eq.  (25). 


The  convergence  of  the  velocity  field  is  assessed  by  testing  whether 


with  analogous  tests  applied  to  the  v  and  w  components  as  well.  €  rv 
and  €av  are  the  relative  and  absolute  convergence  criteria  applied  to  the 
velocity  field;  typical  values  used  are  €  rv  =  10"^  and  €av  =  10"^.  If  Eq. 
(30)  is  not  satisfied,  return  to  Step  3  above  for  another  iteration  cycle. 

Once  Eq.  (30)  is  satisfied,  the  (n+l)st  iterate  is  accepted  as  the  new 
solution  at  (i+1).  The  variable  arrays  are  updated  accordingly,  and  die 
divergence  at  the  center  of  each  cell  is  calculated  and  stored.  If  desired, 
printout  of  selected  portions  of  die  solution  is  made  at  this  stage. 

If  a  predetermined  stop  time  or  a  steady-state  condition  ha:,  been  reached, 
the  calculation  is  halted.  Otherwise,  control  returns  to  Step  2  for 
advancement  through  the  next  time  cycle. 


3.3  Results 


L/H  *  0.3, 
the  way  to 
32  grid  cel 
stretching  ] 
grid  lines  a 
the  endwal 
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conditions  chosen  for  initiai  study  were  H/R  *  3.0,  Re  *  10^,  F  *  0.6, 
and  pps  a  0.  For  these  conditions  the  flow  remains  in  Stage  1  (Figure  1)  all 
solid-body  rotation,  i.e.,  neither  endwali  is  intersected.  N3  =  22  and  NK  a 
11s  were  used  in  the  radial  and  axial  directions,  respectively,  and  the  grid. 
|>arameters  in  Eqs.  (9)  and  (10)  were  b  a  c  »  1.032.  These  choices  put  the 
long  k  a  2  and  j.  *  NJ-l  at  a  distance  6.003  x  10"^  and  2.9  x  10"^  away  from 
1  and  sidewall,  respectively.  The  parameters  oc  and  fj  were  both  zero, 
ng  to  centered  differences. 


Usinj  various  values  of  at,  we  unfortunately  have  not  been  able  to  run  such  a 
case  to  completion.  Typically  the  calculation  proceeds  satisfactorily  for  a  short  time, 
but  with  a  progressively  slower  convergence  of  the  inner  relaxation  for  the  pressure 
field;  the  latter  eventually  fails  to  converge  by  the  maximum  allowed  iteration  cycle 
(ordinarily  >0).  Since  the  pressure  field  is  determined  in  such  a  way  as  to  enforce  the 
continuity  equation,  diagnostics  were  added  to  the  program  to  test  whether  Eq.  (20)  is 


f 


being  accurately  satisfied.  •  These  include  evaluation  of  the  largest  absolute  value  of 
the  divergence  over  the  whole  flow,  where  it  occurred,  and  the  root  mean  square  value 
over  all  the  grid  cells.  Also  evaluated  is  the  overall  mass  flux  integral,  which  in 
dimensionless  form  can  be  approximated  by 


•  J (V  •  V" )  d  Y  »  2  7T  J y  (v-v  )  rdirdi- 

"X  * 

=  2  7r  J J'cv-ir)  r  d  f 


or  2 


Satisfaction  of  Eq.  (20)  thus  assures  local  continuity  while  the  magnitude  of  m  is  a 
measure  of  global  volume  conservation. 

With  a  relatively  large  time  step  of  0.5,  when  convergence  breaks  down  the 
maximum  absolute  divergence  and  its  ron.s.  value  are  approximately  S  x  10"*  and  6  x 
lO"3,  respectively,  and  m  is  about  >5  x  10*4.  These  figures  are  unacceptably  large; 
for  comparison,  the  original  MAC  code  (Ref.  16,  p.  86)  required  that  the  maximum 
divergence  should  not  exceed  3.5  x  10"3.  With  a  much  smaller  time  step  of  0.01,  the 
maximum  and  r.m.s.  values  when  convergence  breaks  down  are  2  x  10" 3  and  3  x  10 
respectively,  with  m  about  1  x  10"^.  While  tnese  values  are  more  acceptable,  the 
troublesome  fact  is  that  they  appear  to  be  steadily  rising,  rather  than  Having  levelled 
off.  Evidently  a  source  distribution  is  somehow  being  numerically  introduced  into  the 
solution.  Presumably  the  convergence  criteria  used  in  Eq.  (29)  will  bear  directly  on 
the  magnitude  of  |  •  But  the  criteria  cited  in  Section  3.7  are  if  anything  even 

more  stringent  than  the  €rp  *  2  x  10~4  of  the  MAC  code  (Ref.  16,  p.  90). 

The  possibility  of  latent  algebraic  or  FORTRAN  errors  is  always  present,  but  at 
this  point  must  be  viewed  as  remote.  More  likely  is  the  possibility  that  the  present 
semi-implicit  method  may  be  susceptible  to  instabilities  not  present  in  an  explicit 
scheme  (Roache^,  p.  195).  A  review  of  the  pertinent  literature  shows  that  Pracht^? 


22.  Pracht,  W.E.,  "  A  Numerical  Method  for  Calculating  Transient  Creep  Flows,"  3. 
of  Computational  Physics,  Vol.  7,  46-60,  1971. 


and  Aziz  and  Heliums  23  were  the  first  to  apply  implicit  methods  to  the  full  Navier- 
Stokes  equations.  Pracht's  interest  was  primarily  in  low  Reynolds  number  flows,  Re  £  1, 
for  which  the  usual  stability  restriction  on  At  for  an  explicit  treatment  of  the  viscous 
terms  is  prohibitively  small.  He  was  able  to  successfully  calculate  such  flows  with  a 
larger  At  by  treating  just  th*  diffusion  and  pressure  gradient  terms  implicitly;  convective 
terms  were  still  evaluated  at  the  previous  time.  Our  interest  here  is  primarily  in 
higher  Re,  and  an  implicit  treatment  of  all  spatial  derivatives.  The  first  such  attempt 
was  that  of  Aziz  and  Heliums 23,  who  applied  an  Alternating-Direction-Implicit  (ADI) 
scheme.  They  found  that  while  this  scheme  worked  well  with  a  stream  function- 
vorticity  formulation  of  the  equations,  when  cast  in  terms  of  primitive  variables  the 
equations  required  prohibitively  small  steps  in  order  to  satisfy  continuity.  They  attributed 
this  to  the  highly  nonlinear  ccjpling  between  the  pressure  and  momentum  equations. 

Thus  for  confined  flows  the  studies  by  Aziz  and  Heliums^  and  Kitchens  ^  indicate 
the  stream  f unction- vorti city  formulation  as  the  method  of  choice,  since  it  allows 
satisfaction  of  the  continuity  equation  ab  initio.  In  the  present  problem  the  presence 
of  the  free  sirface  does  not  easily  admit  such  an  approach,  however.  Further,  hindsight 
suggests  that  the  primitive  variable  approach  taken  in  Ref.  23  may  have  been  flawed 
in  two  respects.  The  first  is  that  a  conventional,  as  opposed  to  a  staggered,  grid  was 
employed;  as  noted  in  Sections  3.3  and  3.5,  this  leads  to  an  inconsistency  between  the 
momentum  and  pressure  equations.  Mare  importantly,  the  dD/dt  term  in  the  pressure 
equation  needed  to  stabilize  the  growth  of  the  divergence  field  was  omitted. 

Mare  recently,  Shadday,  et  al^4  have  studied  the  steady  flow  in  a  rapidly  rotating 
cylinder  with  a  differentially  rotating  endcap  as  the  asymptotic  solution  to  the  full 
transient  equations.  First  order  time  differencing  was  used.  To  eliminate  an  otherwise 
severe  restriction  on  At,  they  found  it  advantageous  to  treat  the  coriolis  and  pressure 
gradient  terms  implicitly.  The  remaining  spatial  derivatives  were  treated  explicitly  on 
a  staggered  grid,  arid  the  term  was  retained  in  the  pressure  equation. 


23.  Aziz,  K.  and  Heliums,  3.D.,  "Numerical  Solution  of  the  Three-Dimensional 
Equations  of  Motion  for  Laminar  Natural  Convection,"  Physics  of  Fluids,  Vol.  10, 
No.  2,  314-324,  1967. 

24.  Shadday,  M.A.,  Ribando,  RJ.  and  Kauzlarich,  3J.,  "Flow  of  an  Incompressible 
Fluid  in  a  Partially  Filled,  Rapidly  Rotating  Cylinder  with  a  Differentially  Rotating 
Endcap,"  3.  of  Fluid  Mechanics,  Vol.  130,  203-218,1983. 


The  only  successful  attempts  at  using  an  Implicit  treatment  for  all  spatial 
derivatives  in  the  primitive  variable  equations  appear  to  be  those  in  Refs.  25-27.  Hodge, 
et  al.25  employed  first-order  accurate  backward  differences  for  the  time  derivative. 
Spatial  derivatives  were  all  approximated  to  second  order,  with  three-point  "upwind" 
differences  used  for  the  convection  terms.  A  point  SOR  method  was  used  to  iterate 
both  the  velocity  and  pressure  equations  at  each  time  step.  Hegna26  later  extended 
this  scheme  to  incompressible  turbulent  flows.  Despite  the  implicit  nature  of  the 
algorithm,  a  time  step  At  <  0.001  was  found  necessary  to  adequately  follow  the 
transient  solution,  resulting  in  CPU  times  of  2-4  hours  (CDC  Cyber  750). 

Moitra27  applied  a  similar  scheme  to  the  three-dimensional  incompressible 
equations.  He  improved  the  time  accuracy  by  employing  second-order  three-point 
backward  time  differencing  in  the  momentum  equations  (as  done  here  in  Section  3.4); 
centered, differences  were  used  for  all  spatial  derivatives,  which  were  treated  implicitly. 
Inexplicably  however,  only  first-order  accurate  differencing  was  applied  to  the  *o/*t 
term  in  the  Poisson  equation  for  the  pressure.  This  produces  an  inconsistency  between 
the  momentum  and  pressure  equations,  in  addition  to  that  introduced  by  his  use  of  a 
conventional,  rather  than  staggered,  grid.  Perhaps  for  these  reasons,  his  solutions  were 
found  susceptible  to  "unusually  high  magnitudes  of  the  divergence."  To  dampen  this 
behavior,  he  introduces  a  spatially  varying  artificial  viscosity  proportional  to  the 
magnitude  of  the  local  divergence.  To  what  degree  this  contaminates  the  numerical 
results  remains  an  open  question. 

References  23-27  did  not  have  a  free  surface,  with  the  exception  of  Ref.  24; 
even  then,  Shadday,  et  al.  assumed  the  free  surface  boundary  conditions  could  be  applied 
along  a  surface  of  constant  radius,  independent  of  time.  Hence  none  of  the  above  had 
the  complication  of  a  moving  surface  boundary  condition  to  contend  with.  For  example, 
this  is  what  prevents  us  from  using  a  more  efficient  direct-solver  or  spectral  method 


25.  Hodge,  3.K.,  Stone,  A.L.  and  Miller,  T.E.,  "Numerical  Solution  for  .  Airfoils  near 
.  Stall  in  Optimized  Boundary-Fitted  Curvilinear  Coordinates,"  AIAA  1*,  Vol.  17, 

No.  5,  458-464,  May  1979. 

26.  Hegna,  H.A.,  "The  Numerical  Solution  of  Incompressible  Turbulent  Flow  over 
Airfoils.  AIAA  Paper  81-0047,  19,th  Aerospace  Sciences  Meeting;  1981. 

27.  Moitra,  A.,  An  Implicit  Solution  of  the  Three-Dimensional  Navier-Stokes  Equations 
for  an  Airfoil  Spanning  a  Wind  Tunnel.  Ph.D.  Thesis,  Dept,  of  Aerospace 
Engineering,  Mississippi  State  University,  1982. 


for  inverting  the  Poisson  equation  for  pressure;  such  methods  are  readily  adaptable  only 
to  flows  with  relatively  simple,  stationary  boundaries. 

The  present  relaxation  procedure  for  the  pressure  attempts  to  enforce  the 
continuity  equation  by  feeding  back  into  the  solution  values  of  the  divergence  at  time 
levels  (i)  and  (i  -  1).  Divergence  values  based  on  the  new  velocities  at  (i  +  1)  do  not 
appear  because  they  are  a  priori  assumed  to  be  zero,  in  the  spirit  of  the  MAC  rpethod 
(cf.  Section  3.5).  Alternative  methods  have  been  proposed  which  use  information  about 
the  divergence  at  time  0  +  1).  For  example,  Hirt  et  al.18  compute  the  divergence  in 
each  cell  based  on  the  latest  velocities.  Then  the  pressure  at  the  cell  center  is  changed 
by  an  amount  proportional  to  minus  the  divergence.  Hence  a  positive  (negative)  value 
of  D  will  produce  a  negative  (positive)  increment  in  p  which  acts  to  drive  D  toward 
zero.  This  recipe  is  followed  from  cell  to  celi  in  the  flow;  a  constant  of  proportionality 
greater  than  one  is  often  introduced  to  overrelax  the '  solution.  Such  a  procedure  can 
be  shown  to  be  analogous  to  solving  a  Poisson  equation  for  the  pressure.  It  has  the 
dual  advantages  over  Eq.  (A-4)  of  having  a  much  simpler  right  hand  side,  and  its  direct 
use  of  the  current  divergence,  as  opposed  to  that  at  previous  times.  Whether  it  could 
successfully  be  incorporated  with  the  present  line  relaxation  for  the  velocity  field  will 
require  further  analysis. 

To  summarize,  we  have  not  yet  been  successful  in  our  application  of  the  PCM  I 
algorithm  to  spin-up  in  a  partially-filled  cylinder.  Our  review  of  the  pertinent  literature 
has  turned  up  several  examples  of  the  successful  application  of  other  implicit  methods 
to  the  incompressible  Navier-Stokes  equations.2^27  The  differencing  scheme  of 
Moitra-27  most  closely  parallels  that  used  here,  although  significant  differences  remain 
in  the  grids  used  and  the  methods  employed  to  solve  the  algebraic  equations. 
Interestingly, '  his  results  also  tended  to  exhibit  unacceptably  large  values  for  the 
divergence  of  the  velocity  field.  We  are  still  in  hopes  that  this  instability  can  be 
removed  without  resorting  to  the  explicit  introduction  of  artificial  viscosity  used  by 
Moitra,  as  the  latter  may  unduly  contaminate  the  results. 

In  any  event,  we  have  found  no  a  priori  reason  to  doubt  the  validity  of  the 
present  approach.  Nevertheless,  the  connected  difficulties  of  slow  convergence  of  the 
pressure  iterations  and  the  inaccurate  satisfaction  of  the  continuity  equation  need  to 
be  addressed.  In  the  future  we  hope  to  study  whether  a  pressure  correction  scheme 
based  on  the  latest  divergence  iterates  could  successfully  overcome  these  problems. 


4.  CONCLUSIONS 


Initially  a  simplified  numerical  model  was  developed  for  the  axisymmetric  spin- 
up  of  fluid  in  a  partially-filled  cylindrical  cavity.  The  analysis  represents  an  extension 
of  the  earlier  treatments  by  Wedemeyer,12  and  Goiler  and  Ranov,H  to  those  cases 
where  the  liquid  free  surface  may  intersect  one  or  both  endwalis.  Earlier  estimates  of 
the  Ekman  layer  pumping  of  the  secondary  flow  are  modified  heuristically  for  situations 
where  the  layer  no  longer  covers  the  entire  wall.  Also,  due  to  the  very  steep  free 
surface  contour  in  the  latter  stages  of  spin-up,  it  was  found  advantageous  to  develop 
the  free  surface  equations  in  an  axial,  rather  than  radial,  coordinate  frame. 

A  computer  program  was  used  to  solve  the  governing  nonlinear  equations  using 
a  straightforward  finite-difference  algorithm.  The  code  predicts  both  the  azimuthal 
velocity  distribution  with  radius,  and  the  free  surface  contour,  as  functions  of  time. 
The  results  exhibit  good  agreement  with  the  time-resolved  experimental  data  for  the 
surface  contour  taken  by  Goiler  and  Ranov.  Their  data  were  taken  only  for  situations 
where  the  free  surface  did  not  touch  the  endwalis.  At  present,  there  are  no  quantitative 
data  against  which  the  simplified  analysis  can  be  compared  for  cases  where  one  or 
both  endwalis  are  intersected. 

Nevertheless,  tbo  following  qualitative  conclusions  have  been  drawn  from  the 
theory.  Plotting  the  fluid  angular  momentum  deficit  vs.  Re"*/Z  t  appears  to  correlate 
the  numerical  data  reasonably  well.  Such  plots  indicate  that  the  angular  momentum 
transfer  follows  a  simple  exponential  behavior  in  time.  For  Re  in  the  range  103  -  103, 
the  growth  rate  appears  uniform  up  to  the  point  where  the  bottom  endwall  is  intersected. 
After  this,  exponential  behavior  is  still  exhibited,  but  at  a  reduced  rate,  reflecting  the 
diminishing  influence  of  the  bottom  Ekman  layer  which  is  primarily  responsible  for  the 
secondary  flow.  The  magnitude  of  the  reduction  in  spin-up  rate  diminishes  with 
decreasing  Reynolds  number. 

To  validate  the  simplified  model  in  a  quantitative  sense,  it  is  desireable  that  its 
predictions  be  compared  with  suitably  designed  laboratory  experiments  and/or  more 
refined  numerical  calculations.  As  an  attempt  to  generate  the  latter,  the  second  phase 
of  this  program  was  devoted  to  developing  a  numerical  model  based  on  the  full, 
axisymmetric  Navier-Stokes  equations.  In  an  effort  to  avoid  the  severe  time  step 


restrictions  characteristic  of  explicit  algorithms,  the  semi-implicit  Predictor-Corrector 
Multi  pie- Iteration  (PCM  I)  finite  difference  scheme  was  applied  to  the  equations  in 
primitive  variable  form.  Unfortunately,  we  have  thus  far  only  been  able  to  march 
the  calculation  for  a  short  time  before  its  convergence  breaks  down.  This  seems  to 
be  related  to  an  unacceptably  large  growth  in  the  divergence  of  the  velocity  field, 
which  should  vanish. 

Our  review  of  the  relevant  literature  has  turned  up  several  examples  of  the 
successful  application  of  other  implicit  methods  to  the  incompressible  Navier-Stokes 
equations.2^-27  Qut  each  differs  in  sufficient  detail  from  the  present  scheme  so  as  not 
to  offer  much  guidance  on  how  to  resolve  the  difficulties.  One  complication  we  face 
which  has  not  been  previously  addressed  by  other  implicit  methods  is  the  presence  of 
a  moving  free  surface.  In  any  event,  we  have  found  no  a  priori  reason  to  doubt  the 
ultimate  applicability  of  the  PCM  I  method  to  this  problem.  Nevertheless,  the  connected 
difficulties  of  slow  convergence  of  the  pressure  iterations  and  the  inaccurate  satisfaction 
of  the  continuity  equation  need  to  be  addressed.  In  the  future  we  hope  to  study  whether 
a  pressure  correction  scheme  based  on  the  latest  divergence  iterates  could  successfully 
overcome  these  problems. 
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grid  stretching  parameters  in  Eqs.  (9)  and  (10) 
divergence  of  the  velocity  field 
Froude  number,  (XLR)2/gH 

axial  acceleration  (  =  gravity  in  a  laboratory  frame) 

height  of  cylinder 

time  index 

radial  grid  index 

axial  grid  index 

initial  fluid  level  before  spin-up 

dimensionless  fluid  angular  momentum  about  the  cylinder's  axis 

number  of  radial,  axial  grid  cells,  respectively 

dimensionless  pressure,  p/p  (aR)2 

dimensionless  free  surface  pressure 

dimensional  cylindrical  coordinates 

radius  of  cylinder 

T/R 


dimensionless  radii  at  which  endwalls  are  intersected  (Figure  1) 
radial  coordinate  description  of  the  free  surface  during  Stage  3 
Reynolds  number,  H.R2/V 
At 


dimensionless  spin-up,  time 

dimensionless  velocities  in  the  (r,  •  ,  z)  directions,  respectively,  normalized 
by  AR 

z/H  in  simplified  model;  z/R  in  Navier-Stokes  code, 
axial  coordinate  description  of  the  fret  surface  during  Stages  1  and  2b 
parameters  used  to  adjust  degree  of  upwind  differencing,  Eqs.  (26),  (A-l)- 
(A-3) 

relative  and  absolute  error  criteria  use|d  in  Eqs.  (29)  and  DO), 
kinematic  viscosity 
density 

axial,  radial  coordinates  in  the  computational  plane 
final  angular  velocity  of  Cylinder,  rad./  sec. 
dimensionless  angular  velocity  of  cylirnter,  normalized  by  -fL 
dimensionless  local  fluid  angular  velocity,  v/r  *  S/a 
indicates  a  dimensional  variable 
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APPENDIX 


We  quote  here  the  full  difference  equations  for  the  velocity  components  and 
pressure,  which  are  somewhat  lengthly  to  include  in  the  main  text.  The  differencing 
procedure  for  the  u,  v,  and  w  equations  is  described  in  Section  3.4,  and  the  results  are: 


u  equation 
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In  writing  Eq.  (A-l),  updated  values  have  been  used  as  soon  as  they  become 
available.  Since  the  grid  is  swept  from  bottom  to  top,  this  means  using  the  (n+1) 
iterate  for  velocities  in  cells  centered  on  row  k  and  below;  above  this  line,  values  from 
the  (n)th  iterate  are  used.  This  scheme  also  saves  considerable  storage,  since  the  (n+1) 
and  (n)th  iterates  need  not  be  stored  simultaneously;  i.e.,  a  single  array  is  used  with 
the  former  replacing  the  latter  as  soon  as  it  is  calculated. 


v  equation 
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The  above  equations  are  tridiagonal  in  the  (n+ 1)  iterate  values  on  each  row  k. 

The  difference  equation  for  the  pressure  is  obtained  'from  a  linear  combination 
of  Eq.  (A-l)  and  (A- 3),  as  described  in  Section  3.3.  The  full  result  is,  for  ec  =  0, 
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